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Abstract 



We derive an analytic expression for the average difference between the forces 



T3 

on the ions in a Car-Parrinello simulation and the forces obtained at the same 

ionic positions when the electrons are at their ground state. We show that for 

common values of the fictitious electron mass, a systematic bias may affect 

the Car-Parrinello forces in systems where the electron-ion coupling is large. 

^^ ■ We show that in the limit where the electronic orbitals are rigidly dragged by 



the ions the difference between the two dynamics amounts to a rescaling of 
the ionic masses, thereby leaving the thermodynamics intact. We study the 



examples of crystalline magnesium oxide and crystalline and molten silicon. 

We find that for crystalline silicon the errors are very small. For crystalline 

MgO the errors are very large but the dynamics can be quite well corrected 

^ ■ within the rigid-ion model. We conclude that it is important to control the 
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effect of the electron mass parameter on the quantities extracted from Car- 
Parrinello simulations. 



I. INTRODUCTION 

Since its introduction in 1985 the Car-Parrinello (CP) methoda has increasingly been used 
to study an ever wider range of problems in the dynamics and thermodynamics of solids 
and liquids under various conditions and in studying the dynamics of chemical reactions. 
CP is an efficient method to solve the Kohn-Sham equations of Density Functional Theorya 
"on the fly", as the electronic ground state evolves due to changing ionic positions. It is 
based on the introduction of an additional inertia associated with the electronic orbitals, 
which are evolved as classical degrees of freedom along the ionic Molecular Dynamics (MD) 
trajectory. Its popularity stems from its efficiency relative to full Born-Oppenheimer (BO) 
dynamical methods, where the electronic orbitals are forced to be in the ground state for 
each ionic configuration, and from the observation that apart from small fluctuations which 
average out on a femtosecond timescale the physical quantities which are extracted from it 
are indistinguishable from BO dynamicsu. 

It is known that the introduction into the electronic system of a fictitious inertia intro- 
duces differences into the dynamics relative to the BO dynamics. However these differences 
have never been fully quantified or analysed in the context of the appropriate theoretical 
framework. In this paper we carefully examine the relationship between CP and BO dy- 
namics and show how the difference between them scales with the value of the fictitious 
mass parameter, /x. We derive analytic expressions for such differences in the limit of the 
fictitious kinetic energy associated with electronic orbitals being a minimum. As has been 
proposed previously^, we show how these errors may be corrected in the ideal case where 
the electron-ion coupling can be modeled as a rigid dragging of localized atomic orbitals and 
show its application to specific examples. We suggest that quantities extracted from CP 
simulations should always be checked against their possible dependence on the value of /i. 

The central issue is the following: the classical motion of the electronic orbitals in CP 
can be thought of as being made up of two components. One component consists of fast 



oscillations with period equal to or faster than r e ~ 27Ta/ /T/ 'E g , E g being the lowest elec- 



tronic excitation energy and // the fictitious electronic massa. The second component is 
the unavoidable but "adiabatic" response of the electronic orbitals to the ionic dynamics, 
whose shortest characteristic time we denote by n (r e <C Tj). It is generally thought that 
keeping the time scales (or frequency spectra) of the two components well separated (i.e. by 
reducing //) ensures a correct adiabatic decoupling and guarantees that the CP dynamics is 
a faithful representation of the BO dynamics. We show here that the adiabatic decoupling 
is a necessary but not sufficient condition for the accuracy of CP dynamics. In particular, 
the fictitious inertia also causes the slow component of the electronic dynamics to exchange 
momentum and energy with the ions, yielding a departure of the CP forces on the ions from 
the BO ones for large values of /i. 

We begin by outlining the relevant theory and deriving an expression for the difference 
between CP forces and BO forces in the limit of minimum fictitious electronic kinetic energy. 
We then consider the simplified model of rigid ions and derive the error in the forces for this 
system. We show how the thermodynamics and the dynamics can be corrected for systems 
in which the rigid-ion model provides a good description of the electronic structure. We 
illustrate these theoretical considerations by studying the examples of Magnesium Oxide 
and Silicon. 

II. THEORY 

The Car-Parrinello method makes use of the following classical lagrangean : 

Lcp = £>ftWi> + \Y, Ml *Z ~ E [ty*h i R /}] (!) 

i I 

to generate trajectories for the ionic and electronic degrees of freedom via the coupled set 
of equations of motion 
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where Mi and R/ are the mass and position respectively of atom J, \ipi) are the Kohn-Sham 
orbitals which are allowed to evolve as classical degrees of freedom with inertial parameters 
fii, and E[{ipi}, {R/}] is the Kohn-Sham energy functional evaluated for the set of ionic 
positions {R/} and the set of orbitals {ipi}. The functional derivative of the Kohn-Sham 
energy in (^) is implicitly restricted to variations of {ipi} that preserve orthonormality. 

We wish to compare the dynamics of ions evolved with this method with the true BO 
dynamics. For this purpose, we decompose the CP orbitals as 

|V>;> = \4 0) ) + Wi) (4) 

where \ip\ ) are the ground state (BO) orbitals which are uniquely defined for given ionic 
coordinates as those that minimize E[{if)i}, {R/}] . This allows us to consider separately the 
evolution of the instantaneous electronic ground state and the deviations of the CP orbitals 
from that ground state. 

A preliminary interesting observation follows from such a decomposition: strictly speak- 
ing, the CP equations do not reduce to the BO equations in the limit of vanishing 8ipi, for 
any finite value of fi. In fact, because of the dependence of the ionic coordinates on time, 
the BO orbitals have a non vanishing "acceleration" given by 

m = TRi^ 1 + YRjR ^f ] (5) 

where we have used the fact that 

It now becomes clear that Eqn. ([|) is not compatible with a vanishing departure of the CP 
orbitals from the BO orbitals, since the right-hand side of ([|) would vanish if {ipi] = {^i K 
while the left-hand side would not, by virtue of (|). So the CP orbitals cannot take their 
ground state values unless fi vanishes too. As a consequence of this, the ionic dynamics is 
affected by a bias proportional to /i and, as we will see, to the strength of the electron-ion 
interaction . 



We now wish to explore the consequences that such a departure from the ground state 
has on the instantaneous CP forces Fqp- We thus calculate how CP forces deviate from the 
BO forces Fbo a t a given point in phase space along the CP trajectory. We may write, for 
the a-th cartesian component of the force on atom J: 
^ a <9£[{R 7 }; frfc}] dE[{R t }; {A}} 
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Substitution of equation |3] yields 
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Using the expansion 
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we can write the error in the CP force as 
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dRJ dRJ 

Having established the connection, to first order in 5ipi, between the CP and the BO forces, 
we assume adiabatic decoupling and look for contributions to this difference that do not 
vanish when averaged over time scales longer than the fictitious dynamics of the electrons 
(r e ) but shorter than the time scales of the ionic dynamics (tj). Only these contributions 
are expected to contribute significantly to the ionic dynamicsu. To this end we rewrite Eqn. 
(|), using (|), as 
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and we show that when eq. (|ll|) is averaged over a time scale shorter than t% but longer 

as 



than r e , then 5ipi vanishes. In order to prove it, we re-express equation 11 
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where 
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and |v?j- ) is an eigenvector of IfW with eigenvalue fej . |x^) and K^ l > , by definition, vary 
on the timescale of ionic motion. We consider timescales much smaller than r» = 2n/uji such 
that |%W),_ftTW anc j j^,^ are approximately constant. In this case a solution of [12] is of the 
form 






M S^)=B^'*-^-l- (15) 



where un = ykj/f^i and B is a complex constant. If we assume that cu is very large 
relative to uj,-, then the average value of ((fj\Sipi(t)) over timescales much greater than 2ir/ujj 
is 



(<Pi(t)\&/„) = — I ( n ( T )\Hi(t))dt 
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this means that for timescales which are intermediate between typical timescales of ionic 



and CP orbital motion \ipi) « |0). 



To summarise : If we consider the dynamics of the electronic orbitals to consist of 
an adiabatic response of the electronic orbitals to the ionic dynamics and an independent 
fast oscillating part then, under the assumption that the timescales of the fast component 
are much shorter than the shortest time period in the ionic system, i.e. assuming adiabatic 
decoupling, the average error in the Car-Parrinello forces is given by (using equations [TO], [TT] 
and 0) 



a*? = 2 v „ iS { v b?^™ + E k , kl a{£\ *m\ + mi „ m (19) 

i Z-^f** \/Lu J dm 8Pl /.^ J k Qm 8RZ8RPJ V l ' V ; 



III. THE RIGID ION APPROXIMATION. 

In order to gain insight into the scale of this problem with the CP forces we consider the 
simple example of rigid ions. We assume that each electron is localised around an ion and 
that there is no distortion of a particular ion's charge distribution as it moves in the field of 
the other ions. We can refer each wavefunction %pi to a particular ion as follows 

^ l (r) = <f) Iv {r-R I ) (20) 

Where the electronic states are labelled by an ion index, /, and the index rj labelling the 
electronic state of the ion . The rigidity of the ionic charge distribution means that 

aMr-R,) = _aMr-R,) and AM* - g0 = Q WJ ^j (21) 

OR i dr dRj 



Equation |19| becomes 
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The second term in equation ^2| vanishes due to symmetry, at least assuming an atomic 
charge density with spherical symmetry. The first term may be written in terms of E k v the 
quantum electronic kinetic energy of an electron in state r\ of atom I as 
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where m e is the (real) mass of an electron. Since the ions are rigid the quantum kinetic 



energy associated with each one is a constant and equation ^2] becomes 

AFf = -AM/.R? (24) 

with 

In this case the ionic positions and velocities are updated during a Car-Parrinello simulation 
according to 

(Mj + AMj)i2? = F£ 0/ (26) 

In other words, for systems where the rigid ion approximation is valid, the CP approximation 
amounts simply to a rescaling of the ionic masses. Since the classical partition function 
depends only on the interaction potential, the thermodynamics of the system as calculated 
with a CP dynamics is identical to the thermodynamics of the BO system. The definition 
of temperature will however be affected, because if the actual ionic dynamics in CP is given 
by (p6|) , then the real temperature at which the system equilibrates, at least in the case of 
a microcanonical dynamics for the ions, is given by 

k * T = ^n ^ {Ml + AM '^<) 2 > ( 2? ) 

I,a 

where (■ ■ •) signifies the average over time and iV is the number of atoms. This differs from 
the standard definition by the addition of a term proportional to AM/. The additional term 
in eq. (|2"T| ) can be readily traced to the additional inertia caused by the rigid dragging of the 
electronic orbitals. In fact, using ([]) and (pT|), we can show that this term coincides, within 
the rigid ion model, with the fictitious electronic kinetic energy, when the contribution from 
the dynamics of the Stpi is negligible i.e. 



T el = J><^> = i^AM 7 (K) 2 ) (28) 

i I, a 

In other words if the electronic orbitals move rigidly with the ions the actual inertia of the 
ions in a CP simulation can be obtained by adding to the "bare" ionic inertia the inertia 
carried by the electronic orbitals. This result has been pointed out previously^ and ionic 
masses are commonly renormalized when dynamical quantities are being investigated. We 
now explore the consequences that such a modification of the ionic inertia has on typical 
observables extracted from CP simulations. First, as already mentioned, the correct defi- 
nition of temperature in a microcanonical CP simulation is given by (|2"T|). Similarly, in a 
simulation where temperature is controlled, e.g. through a Nose thermostattU, the quantity 
to be monitored corresponds to the instantaneous value of (|27|). Dynamical observables will 
also be affected by the additional inertia, as already noted in the case of phonons extracted 
from CP-MD in carbon systemsQ'a. In the case of homogeneous systems (a single atomic 
species in which all the atoms are in similar local chemical environments) all dynamical 
quantities can be simply rescaled using the mass correction of eqn. (^) . However, for het- 
erogeneous systems the correction is not always trivial, as different mass corrections apply 
to different atomic species due to different atomic kinetic energies. In practice, we found 
that a convenient and more general way to express the mass correction of ion I is given by 

9<m Jptotal 

AM - - *t < 29 > 

where fi is a dimensionless constant which takes into account the relative contribution of 
species / to the total quantum kinetic energy El° tal . 

IV. SIMULATIONS 

In order to gain more insight into the theory put forward in the previous sections, we 
have performed CP simulations on pressurised magnesium oxide and on silicon. Among the 
insulators (we restrict our analysis to insulators as adiabatic decoupling is less obvious in 
metallic systems and this would complicate considerably our analysis), MgO and Silicon are 
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extremal cases: MgO is a highly ionic system with large quantum kinetic energy associated 
with the strongly localised charge distribution; Silicon on the other hand is a covalent system 
where electron states are much more delocalised. Within our pseudopotential description of 
MgOB, the Is, 2s and 2p states are frozen into the core of Mg whereas only the Is states are 
frozen into the core of O. Since there is very nearly complete transfer of the two 3s electrons 
from Mg to O (inspection of charge density contour plots reveal no evidence of any valence 
charge anywhere except surrounding O sites) the electron quantum kinetic energy may to a 
first approximation be attributed to electronic states localised on oxygen ions. This makes 
MgO an ideal system to study within the rigid ion model since only the oxygen mass will 
be rescaled. As mentioned in the previous section, additional problems arise if one deals 
with more than one electron-carrying species as the quantum kinetic energy must be divided 
between these species. The large quantum kinetic energy of MgO means that the error in 
the CP forces should be large relative to many materials. The simulations of MgO were 
performed at a high pressure (~ 900 kbar) as this enhanced its ionicity. 

Silicon on the other hand is a covalent/metallic system with relatively low quantum 
kinetic energy. As such it should be one of the systems most favourable to the Car-Parrinello 
approximation but least favourable to description in terms of rigid ions. 

A. Technical Details 

We have performed ten different simulations. The technical details are summarised in 
table | 

All simulations were performed with a cubic simulation cell of side L (see table |) under 
periodic boundary conditions and with 64 atoms in the unit cell. We used a plane wave basis 
set with an energy cut off for the wavefunctions of E cut . The Brillouin zone was sampled 
using only the T— point. In each simulation we have used the mass preconditioning scheme 
of Tassone et alB and the parameters jxq and E p in table | are defined as in Ref . [5] . 

Liquid silicon is metallic and so, as suggested by Blochl and Parrinello Q, two Nose 
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thermostats were used to counteract the effects of energy transfer between the ions and 
the \Stpi) due to overlap of their frequency spectra. The values of the parameters used were 
Q e = 21.3 &.u. / a.tom,Ekin.o = 1-65 x 10 -4 a. u. /atom and Qr = 244400 a.u.. These were 
chosen for compatibility with those of Ref . [10] by taking into account the slight increase in 
temperature and scaling accordingly. 

In all simulations, with the exception of simulation 2, the system was first allowed to 
evolve for at least 1 ps and this trajectory was discarded. For simulation 2, this initial 
equilibration time was 0.5 ps. All results reported are taken from the continuations of these 
equilibration trajectories. 

In all simulations, the total quantum kinetic energy of the system, and hence the average 
mass correction fl29|) , varied during the simulation by less than 0.3%. It was therefore taken 
as a constant in further analysis. 

The total energy of all the degrees of freedom (including the thermostats in simulation 
3) was conserved in all simulations at least to within one part in 10 5 . 

B. Results 

In order to check the predictions of the theory developed in Sections II and III, we have 
taken segments of CP trajectories and calculated the true BO forces along these segments 
by putting the electronic orbitals to their ground state with a steepest descent method. We 
look at the instantaneous error in the a th cartesian component of the CP force on atom / 
relative to the root-mean-squared (r.m.s) BO force component, i.e : 

8Ff{t) = f AF " ( ^ : (30) 



3NN C 

and the instantaneous relative error minus the relative error predicted by the rigid-ion model 



[Hf wu = ^m±^m (31) 

3NN~ C 
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where N c is the number of ionic configurations at which the error in the CP forces was 
calculated and Y2 C * s ^he sum over an suc h configurations. The value of AM/ in fl3~T|) is 
determined using the rigid- ion-model expression (|29|), and El° tal was given its average value 
during the simulation. The scaling parameters which were found to give best results for 
silicon and oxygen were fsi = 1.0 and fo = 1.92 respectively. 

We also look at (5Ff) and {SF") corr the r.m.s values of [5Ff] and [SF^] corr over all the 
ions, cartesian components and configurations tested. 

Since the CP forces are affected by a "fast" component whose effect on the ionic dynamics 
is believed to average out on the time scale of the ionic motion, we introduce the quantity 
T(t), defined as 



v(r)= I — y 



fto+r AFrm 



dr (32) 



j;;; +T \AF?(t)\dt 

If we begin our comparison between CP and BO forces at some instant to along the trajectory 
then inspection of r(r) gives a feeling for how large the fast component is. If the errors in 
all the forces of the system oscillate rapidly with an average of zero then T (r) decreases very 
quickly from the value of one at r = to zero at r ~ r e . For systematic errors T(t) should 
decrease gradually from one to zero on a timescale of the order of the period of r». In realistic 
cases r(r) drops from one and levels off to a smaller value for t ~ r e , and then decreases 
gradually to zero for r exceeding r». The value of T(r) on the plateau between r e and t% 
provides a measure of how much the errors calculated in (|30|) and (|31|) are attributable to a 
systematic (i.e. "slow") departure from the BO surface. 

We begin by looking at the forces in Silicon in both the solid at 330 K and the liquid 
at 2000 K (simulations 1,2 and 3). Simulation 1 was preceded by a short run where the 
temperature was set to about 1000 K. Electrons were then relaxed in their ground state 
and the ionic velocities set to zero. This allows the electrons to smoothly accelerate with 
the ions. A microcanonical simulation followed where the ionic temperature reached, after 
a short equilibration, the value of 330 K. This procedure was followed in all the simulations 
reported here, except where discussed. 

13 



In solid Si at 330 K (Fig. [T] ) we find that the standard deviation of the error in forces 
is 0.94%. However, most of this error can be attributed to a rigid dragging of the Si atomic 
orbitals. The standard deviation of the error is in fact reduced to 0.24% after the rigid-ion 
correction ( p6| ) is subtracted. The ~ 30% drop of T(t) (corrected) shown in Fig. |l]b indicates 
that ~ 30% of the residual 0.24% error can be attributed to "fast" oscillations, so that 
the overall average error introduced by the CP approximation, once corrected for the rigid 
dragging and under the assumption that the fast component is not relevant, is less than 
0.2%. 

As has been pointed out previously by Remler and MaddenE3, it is important to begin 
the dynamics with electrons and ions moving in a consistent way as we have done here 
in all simulations except the one we now discuss (simulation 2) and in the case of liquid Si 
(simulation 3). We found that the error in the forces increases substantially if the simulation 
is not started from zero ionic velocities, a procedure that would otherwise have the advantage 
of shortening considerably the time needed to reach thermal equilibrium. Simulation 2 
started from the end of simulation 1, but electrons have been put in the ground state before 
restarting (ionic velocities and positions were instead kept unchanged). Forces were tested 
after 0.5ps from the electron quenching. 

The standard deviation of the error in forces is now 5.7% and the error in the forces as 
corrected according to the rigid ion approximation at 5.68%, is not significantly improved. 
However clearly from inspection of T(t) in Fig. Qd and the sample force component in Fig. 
0c most of this error can be attributed to the high frequency oscillations of the electronic 
orbitals. If we assume that these oscillations do not influence the ionic dynamics, the error 
reduces to about 1.2% for the uncorrected forces and to less than 0.5% for the corrected 
forces. The amplitudes of these oscillations are nevertheless significant and may affect the 
thermodynamics in a way that is not easy to predict. These oscillations clearly originate 
from the initial jerk experienced by the the electrons in their ground state and survive for a 
long time due to the adiabatic decoupling. 

In the liquid the situation is considerably worse than in the crystal. The standard devi- 
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ation of the error is 3.4% which improves only to 3.1% with the rigid-ion correction. There 
do not seem to be high frequency, high amplitude oscillations here despite the simulation 
being started with finite ionic velocities. However, there are oscillations of a lower frequency 
(although still quite high relative to ionic timescales) which are probably due to the presence 
of the Nose thermostat. It may be that the Nose thermostat has the effect of damping out 
the kinds of oscillations seen in Fig. § but the presence of these other oscillations is hardly an 
improvement. This highlights the need for careful choice of parameters for the Nose thermo- 
stat, particularly the value of -EW n ,o- It is not clear how one should choose this parameter in 
general. For example, here we have used a value of Ef-i n ,o compatible with Ref.[10], however 
we note that this is considerably smaller than the value recommended in Ref . [6] which was 
obtained according to a rigid ion model. We have also done simulations using higher values 
of Eki n fi and in all cases the errors in the forces have been greater. It is likely therefore 
that by decreasing further Eki n ,o we might improve further the forces however this has not 
been attempted here. The issue of thermostatting a system in a way which minimizes the 
kinds of errors seen here while accounting for the evolution of the electronic ground state in 
a more general way than is allowed by the rigid-ion approximation has recently been tackled 
by BldchlH. 

We now look at the forces in crystalline MgO with /j,q = 400a. u (Fig. |]). The relatively 
high quantum kinetic energy associated with states attached to the O ions means that, 
according to eq. (p3[), the errors in the forces are considerably larger for the O ions than 
we have seen for Si. The errors in the CP forces have in fact a standard deviation as large 
as 32% . However, when this is corrected as in (p6|) by attributing all the quantum kinetic 
energy to states rigidly following the O ions the standard deviation of the error reduces to 
4.8%. Furthermore, the corrected value of T(t) indicates that about 80% of the error on the 
O forces cancels out after account is taken for the high frequency oscillations, suggesting that 
a more appropriate estimate of the error is ~ 0.4%. The amplitude of the fast oscillations 
is a cause for concern however and since the simulation was begun at zero ionic velocity it 
is not clear how it may be reduced further. 
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1. Temperature 

We focus here only on MgO, as the effects of the electronic dragging are enhanced. 
According to the results of Section III, we should expect a difference between the naive 
definition of temperature and the definition corrected by the electronic dragging, eqn. (p7|). 
In the case of MgO, as noted in the previous section, this correction affects only the oxygen 
atoms, as only a minor amount of electronic charge is carried by the Mg 2+ ion. In Fig. ^ we 
show the behavior of the instantaneous values of the naive and corrected temperatures. The 
corrected temperature exceeds the naive definition by about 500 K. More interestingly, we 
report in Fig. || also the contributions to the temperature of the two atomic species. It is clear 
that the naive definition would imply that the two species are not at thermal equilibrium. 
On the other hand, use of the corrected definition for the oxygen temperature brings the 
temperature of the two species in much better agreement, supporting the conclusion, based 
on the rigid ion model, that thermodynamics can be restored by a simple rescaling of the 



oxygen mass. The mass rescaling, as calculated with (|22|) amounts to AMq ~ 7.5 atomic 
mass units {Mo = 16 atomic mass units). We also report in Fig. |^ the instantaneous 
value of the fictitious electronic kinetic energy, the l.h.s. of eq. (^), and the difference 
between this quantity and the r.h.s. of eq. (p8|), which represents the contribution due to 
the rigid dragging of the electronic orbitals. The difference is very small, implying that 
residual contributions due, for example, to the fast electronic oscillations are negligible in 
MgO compared to the slow dragging of the orbitals. 

2. Phonon Spectra 

We have calculated the phonon densities of states of crystalline Si and MgO by fourier 
transforming the velocity autocorrelation function. In all cases, the first picosecond of the 
simulation was discarded and results obtained by averaging over at least one subsequent 
picosecond. For Silicon the velocity autocorrelation function was calculated on a time domain 



16 



of length 1.2ps and for MgO on a time domain of length 0.5ps. 

In silicon (Fig. |6|) the difference is reasonably small. According to the rigid ion model 
the frequencies should be corrected using 

^corrected = ^CpV 1 + &M / M (33) 

where ujqp is the frequency as extracted directly from the CP simulation. We find that 
for silicon this overestimates by about a factor of two the amount of the correction. This 
small discrepancy may be due to the length of simulation used for calculating the frequency 
spectra or due to a breakdown of the rigid-ion description when jjlq = 800a. u.. It may also 
be that neglecting the effect of the fast oscillations is not be completely appropriate when 
the dragging contribution is small. 

In MgO, as expected the difference is much larger. We calculate the phonon spectra for 
Ho = 400a. u. and for /i = lOOa.u and find large differences between them (see Fig. [?D 
highlighting again how the dynamics depends on the value of \x. The fact that two species 
are involved complicates matters as the mass correction is different for the two species 
(it actually vanishes for Mg). Therefore we should not expect simply a rigid shift of the 
frequencies. However, if the rigid-ion approximation is valid, one may conceive to rescale 
the oxygen mass a priori in eq. (fj) as Mo = Mo — AMo, so that the actual CP dynamics 
expressed in terms of the BO forces, eq. fl26|) , becomes identical to the BO dynamics if 
the rigid ion approximation holds. We have done this for MgO, again for fio — 400a.u. and 
Ho = lOOa.u and we see that the results are much improved. There are only small differences 
in the positions of the peaks and the overall shapes of the curves are very similar. 

3. Dependence of error on \i 

We now try to address the question of how the error in the CP forces depends on the 
fictitious electron mass. According to Eqn. |1^ the error should scale approximately linearly 
with the mass. However this is based on the simplifying assumptions that the oscillations 
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in \5ipi) have a small amplitude and that the \5ipi) do not on average exchange energy with 
the ions, i.e. full adiabatic decoupling is achieved. Fig. |8] shows (5F") and (SF^) corr for 
the oxygen ion for three different values of /zo where (5Ff) corr has been scaled to eliminate 
the contribution of errors from high frequency oscillations by inspection of T(t). 

Since the uncorrected error is dominated by the effect of the displacement of the equilib- 
rium positions of the orbitals from the ground state, this scales approximately linearly with 
fio- The small error which remains after the rigid ion correction has been applied could have 
contributions from many different sources including deviations from the rigid ion description 
and interactions between the ions and the \5ipi). It is also of the order of the fluctuations in 
j^totai c [ ur i ri g the simulation. 

V. DISCUSSION 

The various cases studied in this work have been explicitly chosen because of their "ex- 
tremal" behavior. The oxygen ions in MgO have a very strong electron-ion coupling as shown 
by the large quantum kinetic energy. However, oxygen seems to be well described within 
the rigid-ion approximation and so the thermodynamics are probably quite close to those of 
a BO system. Crystalline silicon is less well described in terms of rigid ions (although still 
remarkably well) but it has a much lower electron-ion coupling so the errors in the forces 
are very small. If we "measure" the departure of the CP dynamics from the BO dynamics 
in terms of AM/M, with AM defined as in (|26|) , then it appears that the elements where 
the departure is expected to be larger are located in the upper right of the periodic table, 
because they combine a low atomic mass with a large binding energy of the valence electrons 
(and thus a large quantum kinetic energy). Transition metals may also be strongly affected, 
because of the large number and strong localisation of the <i-electrons. However, the higher 
the localisation of the orbitals, the higher the chances that the description of the electronic 
dynamics in terms of rigid orbitals is correct. The large departure observed in the case of 
MgO suggests that a proper assessment of how much the CP forces differ from the BO ones 



is mandatory in most systems. This can be achieved by either calculating the BO forces for 
selected ionic configurations, or by performing simulations for different (smaller) values of 
fi, and checking how the results scale with decreasing jjl. If the departure is large then it is 
likely that in many cases the CP forces can be brought into good agreement with the BO 
ones by simply rescaling the ionic masses. 

Additional complications may arise when the dynamics lead to fundamental changes in 
the electronic structure. The first and second derivatives of the electronic orbitals with 



respect to the positions of the ions, which appear in equation [19], may become relevant in 
regions of phase space where the electrons play a significant role. For example, if charge 
transfer between ions occurs, or if a substantial rearrangement of the electronic orbitals 
takes place, as in a chemical reaction, then the simple method of rescaling the ionic masses 
will no longer work. 

If one is to judge the quality of a CP simulation by the errors in the forces as we have 
largely done here, then a question which needs to be addressed is to what extent these errors 
manifest themselves as errors in the properties of interest in the simulation. It is likely that 
random high frequency oscillations in the forces such as those due to the dynamics of the 
5ipi have no discernible effect on the thermodynamics of the system if such oscillations are 
small. The magnitude of the oscillations seen here in the case of MgO is a cause for concern 
however. 

A further problem may arise at high temperatures where, if there is electron-ion coupling, 
energy can pass between ions and electronic orbitals and the ionic and orbital subsystems 
may begin to equilibrate leading to an increase in the fictitious electronic kinetic energy. 



This can be seen by inspection of equation 12, where the ionic terms may vary over short 



timescales if the ionic kinetic energy is very high and the fictitious mass /i too large. It has 
previously been supposed^ that such an irreversible increase of kinetic energy was due to 
overlap of the ionic frequency spectrum with that of the electrons due to a small or vanishing 
bandgap, E g . However a small bandgap is a sufficient but not a necessary condition for this 
energy transfer to occur. 
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All of the effects discussed in this paper are dependent on the choice of the fictitious mass 
parameter, //, and by reducing this parameter all thermodynamic and dynamic properties 
of a simulation may be brought arbitrarily close to those in a Born-Oppenheimer system. 
However, a reduction of [i has the drawback that the time step required to integrate the 
equation of motion for the electronic orbitals is reduced thereby decreasing the computational 
efficiency of the method. By checking how the property of interest in a simulation scales 
with \x one can control the level of approximation with which it is calculated. 

VI. CONCLUSIONS 

Under the assumption that high frequency electronic oscillations (i.e. the dynamics of 
the \Sipi)) are small and independent of ionic motion, we have shown that Car-Parrinello 
simulations amount to solving the equation of motion for the ions 



Ml R« = F« BO + 2Y^\YR ff ' 9 M? +y^Rl d f j ^ X I (34) 



We have compared the forces in simulations of Si and MgO for a number of different 
values of \x to the BO forces and found that in the case of Si the errors are small and 
change very slightly the phonon spectrum of the crystal. In MgO we observe very large 
systematic errors in the forces which are however mostly attributable to a rescaling of the 
mass of the oxygen ion thereby preserving the thermodynamics. When corrected for this 
effect the errors are slightly higher than those in crystalline Si but still quite small. The 
phonon densities of states further confirm both the inadequacy of CP at these values of n 
without the rigid-ion correction to describe dynamics and the ability of the rigid-ion model 
to correct the dynamics in MgO. 

We have demonstrated the necessity for checking the dependence of results of CP simu- 
lations on the value of the fictitious mass parameter //. 
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TABLES 





TABLE I. 


Technical Details of the Simulations 






Simulation # 


System 


Temperature 


Po 


E p 


Ecut 


At 


EiW(^IV'i) 


L 


Kelvin 


a.u. 


Ryd. 


Ryd. 


a.u. 


a.u.xlO 4 


a.u. 


1 


Si 


330 


270 


1.0 


12.0 


5.0 


4.36 


20.42 


2 


Si 


330 


270 


1.0 


12.0 


5.0 


4.36 


20.42 


3 


Si (liquid) 


2000 


270 


1.0 


12.0 


10.0 


4.35 


19.8 


4 


MgO 


2800 


400 


2.7 


90.0 


8.0 


66.3 


14.5 


5 


Si 


330 


200 


1.0 


12.0 


5.0 


3.23 


20.42 


6 


Si 


330 


800 


1.0 


12.0 


10.0 


12.92 


20.42 


7 


MgO (M rescaled) 


2800 


100 


2.7 


90.0 


4.0 


16.55 


14.5 


8 


MgO (M rescaled) 


2800 


400 


2.7 


90.0 


8.0 


66.4 


14.5 


9 


MgO 


2800 


200 


2.7 


90.0 


5.65 


33.1 


14.5 


10 


MgO 


2800 


100 


2.7 


90.0 


4.0 


16.55 


14.5 
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FIGURES 




0.02 0.03 0.04 0.05 0.06 
Time [picoseconds] 
FIG. 1. Simulation 1. (a) Distribution among all atoms / and all cartesian components a 

of the percentage errors in the CP forces relative to the BO forces at the same ionic positions, 

100 x 5Ff(t) (full line) and these errors when the forces have been partially corrected according 

to a rigid-ion model, 100 x [5Ff(t)] corr (dashed line), (b) T(t) as defined by Eqn. |3^ for the full 

error in the forces and those as partially corrected according to the rigid-ion model (c) F^q ,F^ p 

and {F(jp — Fg ) (multiplied by a factor of 20 for visibility) for a typical force component. Dots 

indicate the points at which the BO force was calculated (every 5 time steps). 



24 




0.02 0.03 0.04 

Time [picoseconds] 
FIG. 2. Simulation 2, Crystalline Si at 330K when the electrons receive a 'kick' at the beginning 

of the simulation. See caption of Fig. 1 for explanation. (Fq P — F^ ) has not been scaled for 

visibility. 
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0.02 0.03 0.04 
Time [picoseconds] 
FIG. 3. Simulation 3, Liquid Si at 2000 K. See caption of Fig. 1 for explanation. (Fq P —F^q ) 



has not been scaled for visibility. 
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-60 -40 -20 20 40 60 

(b) 1 -, v_^ Relative Error [%] 




0.003 0.006 0.009 

Time [picoseconds] 
FIG. 4. Simulation 4. Forces on the oxygen ions in crystalline MgO at 2800 K. (a) 

and (b) are as in Fig. 1 (c) From top to bottom (Fq P — Fg ),—AMoRf (dotted line), 

(Fq P — Fg + AMoRf),F^ p ,Fg for a typical force component. Dots indicate the points 

at which the BO force was calculated (every time step). 
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1.1 1.15 1.2 1.25 

Time after beginning of simulation [picoseconds] 
FIG. 5. Simulation 1. Dashed lines from top to bottom are : Mg temperature, O temperature, 

T e \ ,and 10 x (T e ; — Tam )- The full line is the Oxygen temperature when it is calculated with a 

mass which is increased by Aq. 



|x =200 a.u. 
|j, =800 a.u. 




400 

Frequency [cm~ : ] 
FIG. 6. Phonon density of states of crystalline Silicon for [Iq = 200 a.u. (simulation 5) and for 

fio = 800 a.u (simulation 6). 
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p. = 100 a.u 




200 



400 



600 
Frequency [cm 



800 

Hi 



1000 



1200 



FIG. 7. Phonon density of states of crystalline MgO for [Xq = 100 a.u. and for hq = 400 a.u. 
with rescaled (simulations 7 and 8) and unrescaled (simulations 10 and 4) oxygen masses. 
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FIG. 8. Scaling of the standard deviation of the errors in the forces on the oxygen ions with 

fj,Q (Simulations 4,9 and 10). {5Ff) corr has been reduced to eliminate cancelling high frequency 

oscillations by inspection of T(t). 
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